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We study the phase structure of QCD at finite temperature and density by numerical 
simulations on a lattice. The most important point for the numerical study at finite density 
is treatment of the sign problem. We propose a method to avoid the sign problem, which is 
based on a cumulant expansion of the complex phase in the density of state method combined 
with the reweighting method. Using the method, we study the critical point terminating a 
first order phase transition line in lattice QCD at high temperature and density. 



It is very important to explore the QCD phase structure to understand the his- 
tory of the universe. We expect that the nature of the chiral phase transition changes 
as the quark number density increases. Moreover, new state of QCD matter may 
appear at high density. However, because the quark determinant is complex at finite 
density, the Monte-Carlo method is not applicable directly for finite density QCD. 
One of the popular methods to avoid this problem is the density of state method. 
We adopt an appropriate physical quantity such as quark number, chiral order pa- 
rameter, gauge action etc., which is denoted by X, and discuss its state density. The 
state density, i.e. the probability distribution function, at finite temperature T and 
quark chemical potential \x q is defined by 



where 5(x) is the delta function, M is the quark matrix, S g is the gauge action, 
and Nf is the number of flavors. The partition function Z is given by Z(T, fj, q ) = 
f W(X,T, fj, q )dX. Once we obtain the probability Ql-ip . expectation values of the 
operator O of X, e.g. (X), ((X 2 — (X)) 2 ), can be evaluated by the following equation; 



Although these equations are rather trivial, the probability distribution of a physical 
quantity, Eq. (|l-ip . is well-defined as a real number even when detM is complex. 

In this report, we discuss the density of state approach combined with the 
reweighting method to investigate the QCD phase structure at high density. In 
the next section, we discuss the state density using the reweighting method. For this 
calculation, we introduce a method to avoid the sign problem in Sec. El Using the 
method, we study the phase diagram. The investigation of the distribution function 
is one of the most primitive approaches to identify the order of phase transitions. We 
expect that two phases coexist at a first order phase transition point. In Sec. HJ we 
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Fig. 1. Left: The plaquette histogram w(P, /3) and the effective potential V e g (P, /3, 0) at [i q = for 
each p. Right: The histograms of complex phase 8 for /iq/T = 1.0 and 2.0 at ft = 3.65.^ 

calculate the distribution function and discuss the order of phase transitions using 
the distribution function. Conclusions are given in Sec. 



§2. Density of state in the reweighting method 



We discuss the density of state method with fixing the plaquette variable (P), 
i.e. lxl Wilson loop, as an example. The distribution function is defined by 
Eq. (jl-ip with X = P. For later discussions, we define the average plaquette P 
as P = —S g /(6/3N s i te ) and the quark matrix M as independent of f3. N s \ te is the 
number of sites, and the parameter j3 = 6/g 2 controls the temperature. When the 
quark determinant is real, the distribution function is given by the histogram of 
P. The plaquette distribution function and V e ff(P) = — InW(P) for p4-improved 
staggered fermions at \i q = obtained in Ref. [1]) are shown in Fig. Q] (left). V e g (P) 
is normalized at the minimum point, and the minimum point moves to right as the 
temperature (T/T c ) or (3 increases. T c is the transition temperature at fi g = 0. We 
denote the distribution function at [x q = as w(P,f3) = W(P, /3,0). 

Because the quark determinant is complex at finite n q , the reweighting method 
is used to obtain the distribution function.^ The partition function is rewritten as 

Z((3, N ) = J W{P,(3, N ) dP = f R(P, N )w{P,P) dP. (2-1) 

Here, R(P, n q ) is the reweighting factor for finite \i q defined by 

UtV - m (detM(n q )) N ! \ 
mpl , $VU5{P' -P)(tetM{^)) N t ^> (detM(O))^ / ( ^ =0) 

H ^,W- fvU6(P'-P)(detM(0)r< (t(P'-P)) { ^ q = 0) [ } 

This R(P, n q ) is independent of j3, and R(P, fi q ) can be measured at any /3. In this 
method, all simulations are performed at fj, q = and the effect of finite \i q is in- 
troduced through the operator det M(fi q )/ detM(0) measured on the configurations 
generated by the simulations at fj, q = 0. 
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Since QCD has the symmetry of charge conjugation, the partition function is 
invariant under a change from \i q to — (j, q , i.e. R(P,—fi q ) = R(P,fi q ). Moreover, 
the quark determinant satisfies detM(— (j, q ) = (det M (fj,*))* . From these equations, 
we get [R(P, fig)]* = R(P,fi*). This indicates that R(P, fj, q ) is real if fi q is real, 
i.e. fi q = n*, and the probability distribution function of the plaquette given by 
W{P,P,ix q ) = R(P,n q )w(P,f3) is real. 

§3. Avoiding the sign problem 

However, a serious sign problem occurs in the calculation of R for large fi q /T. 
The histogram of the complex phase 9 are shown in Fig. [TJ (right) obtained in a 
simulation by p4-improved staggered fermionsP 1 The complex phase of the quark 
determinant is defined by a Taylor expansion; 



0(ji) = iV f Im[lndetM(^) 
1 



n (2n + l)! 

n=0 



Im 



d 2n+1 (mdetM(/i g )) 
0(/z 9 /T)2"+i 



(fJ-q/ T ) 2n+1 - (3-1) 

(m«=o) 



We note that lndetM(// 3 ) is not uniquely defined for complex det M(fj, q ). The 9 
defined in Eq. f)3- 1 j) is not restricted to be in the range —it to it, and the maximum 
value of \9\ is infinite in the large volume limit. Of course, we can restrict the range 
of 9 from — n to ir subtracting 2nn, where n is an integer, in the definition of 9. 
However, this ambiguity does not affect the calculation of (e ). 

In Fig. [TJ (right), the width of the distribution becomes wider as n q /T increases, 
corresponding to the phase fluctuation larger. The expectation value of (e %e \ det M |) 
decreases as the fluctuation of 9 increases, and the expectation value becomes smaller 
than the statistical error when the complex phase fluctuation of the quark determi- 
nant becomes larger than 0(tt) in the Monte-Carlo steps. This is the sign problem 
in the calculation of the reweighting factor. 

To avoid the sign problem, we perform the 9 integration before the integration of 
| det M(Atg)/ det ^"(0)1^ = F in the calculation ofEq. $T2§; (e i0 F) = f(e ie ) F FdF, 
where means the expectation value with fixed F. We then consider the 

following cumulant expansion; 



(e ld ) F = exp 



»Wc--2 ^r + — + —\ + 



(3-2) 



where (9 n ) c is the n th order cumulant, e.g. (9 2 ) c = (6 2 ) F , {^) c = (^)f ~ 
3(9 2 f F , (9 & ) c = (9 & ) F -lh(9 A ) F {9 2 ) F + m{9 2 f F . Note that (9 n ) c = for odd 
n due to the symmetry under 9^—9. Because only the odd-order cumulants are 
the source of the complex phase in (exp(i#))i?, the value of (exp(i#))j? is guaranteed 
to be real and positive from this symmetry if the cumulant expansion converges. 
Although the identity (|3-2p is exact if we consider infinite terms of the expansion, 
there is no source of the sign problem once we eliminate the odd terms. 

As shown in Fig. [TJ (right), the distribution of 9 is well-approximated by a Gaus- 
sian function (dashed line). When the distribution of 9 is Gaussian, the 0(9 n ) terms 
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vanish for n > 2 in Eq. (|3-2p . Hence, the approximation that the higher order cu- 
mulants are neglected except for the first nonzero term is equivalent to the Gaussian 
approximation for the 9 distribution. When one wants to improve the Gaussian 
approximation, it is achieved by adding higher order terms. 

Moreover, the cumulant expansion can be regarded as a power expansion in 
terms of fi q because 9 ~ 0(fj, q ). Therefore, if we take into account the cumulants 
up to the n th order, the truncation error does not affect the Taylor expansion up 
to 0(/i"). The Gaussian approximation corresponds to the leading non-trivial order 
approximation of the Taylor expansion in [x q . 

On the other hand, a careful discussion about the infinite volume (V) limit is 
required.® Because the operator 9 is roughly proportional to V, the n th order cumu- 
lant {9 n ) c may increase as 0(V n ) naively. If this is the case, the cumulant expansion 
does not converge at large V. However, the following argument suggests that the 
convergence property of the cumulant expansion is independent of the volume when 
the correlation length of the system is finite. Note that, since no critical point is 
expected to exist in two-flavor QCD at m q > and fi q = 0, the correlation length 
between quarks is finite. The expansion coefficients of 9 in Eq. (|3-ip are given by 
combinations of traces of products of M _1 , d n M/d(fi q /T) n and so on. For example, 
the first coefficient is given by the trace of N f [M- 1 {dM/d{iJ Jq /T))} and the diagonal 
element of this matrix is the local quark number density operator (~ ^^^(x)) at 
H q = 0. When the correlation length of the local number density operator is much 
shorter than the system size, we may decompose the first derivative term into in- 
dependent contributions from spatially separated regions. The same discussion is 
applicable to higher order coefficients too. 

In this case, one can write the phase as 9 = ^2 X 9 X , where 9 X is the contribution 
from a spatial region labeled by x and these contributions are independent. The 
average of exp(i0) is thus 

x \ x n / 

This equation suggests that all cumulants (9 n ) c m Y^ x {&x) c increase in proportion to 
the volume as the volume increases. Therefore, while the width of the distribution, 
i.e. the phase fluctuation, increases in proportion to the volume, the ratios of the 
cumulants are independent of the volume. The higher order terms in the cumulant 
expansion are well under control in the large volume limit. 

In addition, the complex phase can be decomposed into independent parts when 
we define 6 as 9 = J M?/T Im[dlndet M/d{fi q /T)] { ^ /T) d(n' q /T), as well as Eq. (13~T|) . 

Because 9 is 0((j, q ) and (9 n ) c is 0(/j, q ), the Gaussian approximation is valid 
at small [i q and the higher order cumulants will become visible at large fj, q . The 
application range of the Gaussian approximation in terms of fi q must be checked for 
each analysis by calculating the ratio of cumulants. However, it is expected from 
the argument of the volume-dependence of the ratios that the application range does 
not change once the system size becomes larger than the correlation length. This 
property will enable us to use large lattices. 
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Fig. 2. Left: The derivative of the effective potential dV e g/dP at /3 = 3.65.^ We measured them at 
the peaks of the plaquette histograms in Fig. [1] (left) and interpolated the data by a cubic spline 
method. Right: The derivative of In Zc vs. quark number density computed with a saddle 
point approximation™ T c is the transition temperature at \x q = 0. 



§4. Distribution function and first order phase transition 

In this section, we discuss the order of the phase transition at finite density^ 1 ® 
using data obtained in simulations with the 2-flavor p4-improved staggered quarksP 
m n « 770MeV. The distribution function is expected to be a double-peaked function 
at a first order transition point, i.e. V e s = — In W is a double- well function. It is easy 
to prove useful properties in the plaquette effective potential; From the definition of 
W and P, W(P,(3,n q ) = e^'^^^W (P, /3o, fj, q ) is satisfied under the parameter 
change from /3n to (5. Then, dV e g/dP at different j3 can be estimated by the equation; 

^?CP,/W = ^(P,/3o,M g )-6(/3-/3o^ite, (4-1) 

and d 2 V e g/dP 2 is independent of [3. Therefore, the shape of dV e ^/dP as a function 
of P does not change with [3 up to a /3-dependent constant. 

Using Eq. (|4T|) . we calculate dV e g/dP at \x q = in a wide range of P. Perform- 
ing simulations at many j3 and finding the peak position of the plaquette distribution, 
at which dV eS /dP(P) = 0, the value of dV eS /dP(P) for (3 is given by 6(/3 -f3 )N site . 
This method is much easier than the estimation from the plaquette histogram be- 
cause the range of plaquette value obtained by a simulation with single /3 is narrow. 
(See Fig. [1] (left).) The finite density effect of V e g is evaluated calculating R(P,fj, q ) 
with the Gaussian approximation. We plot dV e g/dP instead of V e g itself in Fig. [2] 
(left) for various n q /T. j3 = 3.65 is adopted for these results, however the f3 value can 
be easily changed by Eq. (|4-ip . If the effective potential V e s is a double- well function 
of P, there exists a region of P where the derivative of dV c g/dP is negative. The left 
panel of Fig. [2] shows that the region of d 2 V eS /dP 2 < exists for (n q /T) 2 > 6. This 
suggests that the phase transition becomes first order at high density. The details 
of this analysis are given in Refs. [2]),|3]). 

Finally, we want to mention the distribution function of the quark numberP 
The probability distribution is, in principle, measurable by event-by-event analysis 
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of heavy- ion collisions. The Gaussian approximation is also useful for the calculation 
of the quark number distribution function. The relation between the grand canonical 
partition function Z and the canonical partition function Zq is given by the following 
Laplace transformation; 

Z(T,n q ) = ]T Z c (T,N)e N ^ T = V I Z C (T, P V)eP v ^' T dp. (4-2) 

N J 

where N is the quark number, V is the volume and p = N/V is the quark number 
density. — \hZq — pV p q /T is regarded as the effective potential V e ${p). 

We compute the derivative of In Zc with respect to p by the saddle point ap- 
proximation using the data obtained in Ref. Q. If the distribution function is a 
double-peaked function, the derivative of InZc is an S-shaped function. Here, we 
denote p* /T = — (l/V)dlnZc/dp, since p*/T = p/T in the thermodynamic limit. 
This calculation suffers from the sign problem. To eliminate the sign problem, the 
approximation discussed in the previous section is used, i.e. the complex phase factor 
e %e is replaced by exp[— (6 2 )/2\. The details are given in Ref. [5]). 

The result of p* q /T is shown in Fig. [2] (right) as a function of p/T 3 for each 
temperature T/T c . The dot-dashed line is the value of the free quark-gluon gas in 
the continuum theory, p/T 3 = Nf[(p q /T) + {l/-n 2 ){p q /T) 3 \. From this figure, we 
find that a qualitative feature of p* q /T changes around T/T c ~ 0.8, i.e. p q /T in- 
creases monotonically as p increases above 0.8, whereas it shows an S-shape below 
0.8. The behavior at low T is a signature of a first order phase transition. Although 
some approximations are used, the critical value of T is roughly consistent with the 
critical point estimated by the plaquette effective potential using the same configu- 
rations, (T/T c , pL q /T) m (0.76,2.5)!^ The difference between these two results may 
be a systematic error. Further studies are necessary to predict the critical point 
quantitatively, but these results are consistent with our qualitative expectation. 

§5. Summary 

We discussed methods to investigate finite density QCD beyond the low density 
region. A method based on the investigation of an effective potential as a function 
of the average plaquette was proposed introducing an approximation to avoid the 
sign problem, and the existence of the critical point at finite density is suggested 
by simulations with improved staggered quarks. Moreover, it was found that inter- 
esting information about the QCD phase structure at finite density is obtained by 
constructing the canonical partition function for each quark number. 
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